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Abstract. The linear instability of thin, vertically-isothermal Keplerian discs, under the influence 
of axial magnetic field is investigated. Solutions of the stability problem are found explicitly by 
asymptotic expansions in the small aspect ratio of the disc. It is shown that the perturbations are 
decoupled into in-plane and vertical modes. Exact expressions for the growth rates as well as the 
number of unstable modes are derived. Those are the discrete counterpart of the continuous infinite 
homogeneous cylinder magnetorotational (MRI) spectrum. In addition, a weakly nonlinear analysis 
of the MRI is performed. It is shown that near the instability threshold the latter is saturated by the 
stable magnetoacoustic modes. 
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INTRODUCTION 

The stability properties of Keplerian disks have been a focus of intensive investigation 
of theoretical astrophysicists over the last decades, pertaining to the problem of angular 
momentum transfer in accretion disks. Traditionally, results of the analytical study of 
magneto-rotational instability (MRI) in an infinitely long cylinder ([1]; [2]) have been 
adopted for thin disks in order to derive criteria for the spectral stability under various 
conditions ([3]). In addition, it has recently been shown that if the thin disk geometry is 
taken into account both the growth rates as well as the number of unstable MRI modes 
are greatly reduced and are decreasing functions of the disk thickness ([4]; [5]; [6]). The 
aim of the current work is twofold: 1. to provide a complete description of the stable as 
well as of the unstable spectrum of thin disks under the influence of an axial magnetic 
field, and 2. to carry out a weakly nonlinear analysis of the MRI in order to investigate 
its saturation mechanism. 



THE PHYSICAL MODEL OF THIN KEPLERIAN DISKS 

The stability of radially as well as axially stratified rotating plasma in thin vertically 
isothermal discs threaded by a magnetic field is considered. Viscosity, electrical resis- 
tivity, and radiation effects are ignored. 



Governing equations 



As a first step, all physical variables are transformed to non-dimensional quantities by 
using the following characteristic values ([6]; [7]): 
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Here £2* = [GM c j r\) l l 2 is the Keplerian angular velocity of the fluid at the characteristic 
radius r* that belongs to the Keplerian portion of the disc; G is the gravitational constant; 
M c is the total mass of the central object; c is the speed of light; 4>* is the characteristic 
value of the gravitational potential; the characteristic mass and number density equal to 
the ion mass and number density, m* = m,- and = n\. The characteristic values of the 
electric current density and electric field, j* and E* have been chosen consistently with 
Maxwell's equations; cs* is the characteristic sound velocity; T* = T(r*) is the charac- 
teristic temperature; the characteristic magnetic field is 5* = 5 z (r*). The dimensional 
equilibrium temperature T(r) and B z (r), are free functions. 

The resulting dimensionless dynamical equations for vertically isothermal discs are: 
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^ + VxE = 0,V-B = 0, (4) 
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E = —V x B, (5) 

P = nT. (6) 

Here VP = c|Vn for vertically isothermal discs, and the dimensionless equilibrium 
sound speed is given by c 2 s = dP/dn = T(r). Standard cylindrical coordinates {r, 6,z} 
are adopted throughout the paper with the associated unit vectors {i r ,ig,i z }; V is the 
plasma velocity; t is time; D/Dt — d /dt + (V ■ V) is the material derivative; 4>(r,z) = 
— (r 2 +z 2 ) -1 / 2 is the gravitational potential due to the central object; B is the magnetic 
field, j = V x B is the current density; E is the electric field; P = P e + Pi is the total 
plasma pressure; P/ = n/7} are the partial species pressures (/ = e,i); T = T e = T t is the 
plasma temperature; subscripts e and i denote electrons and ions, respectively. Note that 
a preferred direction is tacitly defined here, namely, the positive direction of the z axis is 
chosen according to positive Keplerian rotation. The dimensionless coefficients Ms and 
/3 are the Mach number and the characteristic plasma beta, respectively: 



Zero conditions at infinity for the in-plane magnetic field and the number density are 
adopted, namely: 

B r = 0, B e = 0, n = for z = ±°°. (8) 

To simplify the further treatment of Maxwell's equations, both the hydro-magnetic basic 
configuration and perturbations are assumed to be axisymmetric. 

A common property of thin Keplerian discs is their highly compressible motion with 
large Mach numbers ([8]). Furthermore, the characteristic effective semi-thickness of 
the equilibrium disc //* = //(r*) (H = H(r) is the local semi-thickness) is defined so 
that the disc aspect ratio £ equals the inverse Mach number: 
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Thus, the thin disc approximation means 
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The smallness of e means that dimensionless axial coordinate is also small, i.e. z/r* ~ 
£ and consequently the following rescaled quantities may be introduced in 

order to further apply the asymptotic expansions in £ [(similar to [6]; [7]; [9]): 

C = ^ e °,/f(r) = ^~£°, (11) 
where H(r) = cs(r)/Cl(r) is the scaled semi-thickness of the disc. 



Equilibrium configurations 

We start by deriving the steady state solution. It is first noted that the asymptotic 
expansion for the time-independent gravitational potential is given by: 

4>(r, = *(r) +£ 2 <Kr, Q, <J>(r) = -I 0(r, £) = ^C^W + C(£ 2 ), r > 1 » £. 

(12) 

Substituting (12) into Eqs. (2)-(6) and setting the partial derivatives with respect to time 
to zero yield to leading order in £ 

Vj = d*(r) c 2 s {r)dn_ 30(r,Q _ (13) 
r dr ' n dC, dC, 



Thus, the velocity as well as the number density are given by: 

K = o(e 2 ), V e = £°V e (r) + 0(£ 2 ), V z = o(e 2 ), n = £°n = e°iV(r) v(tj), (14) 
where o(e) e, 0(e) ~ e, and 

Vfl(r) = rQ(r), Q(r) = r^'\ vfa) = exp(-r] 2 /2), T] = C/tf (r). (15) 

There are two solutions of eqs. (13) - (15) that describe two quite different equilib- 
rium magnetic configurations, namely one with comparable magnitudes of the axial and 
toroidal components of the magnetic field, and the other with dominant toroidal compo- 
nent. The two equilibria are distinguished by different scaling of the physical variables 
with e. A detailed description of those two equilibria may be found in [6] . Here however, 
we focus on the equilibrium that is characterized by comparable magnitudes of the axial 
and toroidal components of the magnetic field, and in particular will consider the case 
of a pure axial equilibrium magnetic field. As is shown in [6], the results may be easily 
extended to the more general case of a comparable toroidal component. All equilibrium 
variables are written in the leading order in £, and depend on the radial variable only. 
The exceptions are the number density and the pressure that depend on the axial coordi- 
nate in a self-similar manner with radius-dependent amplitude. The axial magnetic field 
as well as the disc thickness and the amplitude factor, N{r), in the number density are 
arbitrary functions of the radial variable. 

To start the equilibrium description it is first assumed that the axial component of 
the magnetic field is of order e°. That assumption together with relations (12)-(15), 
determine the order in e of the rest of the physical variables. The result is given by: 

B r = o(e 2 ), B e = e l B e (r), B z * e°B z (r), (16) 

j r -o{£ ), je = £ Jev) = -e — , j z = e j z = e — • (1 /) 

dr r dr 

It is noted finally that to lowest order in £ the magnetic field configurations under 
consideration do not influence the steady-state properties of the disk. As will be seen in 
the following sections, this situation changes dramatically when small perturbations are 
considered. 



Perturbed thin discs 

In general for the unsteady nonlinear case the dependent variables are scaled in e in 
the following way: 

f(r,t,t) = £ S f(r,0+£ S 'f(r,t,t). (18) 

Here f(r, stands for any dependent variable, the bar and the prime denote equilib- 
rium and perturbed variables; each perturbed variable is characterized by some power 
S', as is summarized in Table 1. 



TABLE 1. The e-scaling of the perturbed variables. 

f = e S f + e s 'f n V r V e V z B r B 6 B z j r jg j\ 
S' 1110 1-1-10 



DYNAMICAL EQUATIONS FOR THE PERTURBED DISCS 

As stated above, a detailed stability study is carried out for zero toroidal magnetic field. 
As shown in [6] the results of such analysis may be extended to the case of comparable 
poloidal and toroidal equilibrium components. 



The reduced nonlinear equations 

The perturbations are subject to the boundary conditions for the in-plane magnetic 
field and number density as follows from (8) 

B' r = 0, B' e = 0, ri = for £ = ±°o. (19) 

Both poloidal and toroidal components of the perturbed magnetic field, B' r and B' z , 
are expressed through the magnetic flux function, , which renders the magnetic field 
divergent free, V • B = 0. 

As the radial coordinate is a mere parameter in the set of the reduced equations, it is 
convenient to replace the physical variables by the following self-similar quantities: 

T = », V = ^ 7 y (20) 



such that the derivatives in the new and old variables are related as follows: 

d d d Id 



(21) 



Finally, a simpler form of the final equations is obtained by introducing the following 
scaled dependent variables: 

v(T,r,r ] ) = ^, v(T,r,r 7 ) = ^ y , b(r,r, rj) = j^r y (22) 

Below for convenience and with no confusion the notation t for the time variable is 
reinstated instead of the new variable T. This yields the following system of equations 
that depend parametrically on the radius: 

1 dv r 1 1 db r dv r 

'2ve--— t^-^— ^— = -v z — , (23) 



&(r) dt p( r ) v{ri) + vdri drj ' 

1 dv e 1 1 1 db e dv e 

+ ^ v r- ir, \ .-./v — ^— = - v z^-, (24) 



d(r) dt 2 r j8(r)v(?]) + V dr] z drj ' 



'<"> 3 .( » )= »*[^ + 4_4 + « ] , (25) 



ft(r) v(77) + v<977 v v(r]) y 2 <9t7 L z /3(r) v(tj) + v J 

i av + 3[Wn)+vK =0 (26) 



1 db r dv r d{y z b r ) 
Cl(r) dt dt] ~~ dt] 



(27) 



1 db e dv e 3 d(v z fr ) 

+ -0 r = =r~~ — , (28) 



Q(r) dt dt] 2 <9r] 

supplemented by the vanishing conditions for the in-plane magnetic field and number 
density at infinity. Here the value of the epicyclical frequency in a Keplerian rotating 
medium, %(r) = Q(r) has been employed; v(t]) is the scaled equilibrium density; j8 (r) 
is the local plasma beta function: 

where the local parameter /3 (r) is proportional to the characteristic plasma beta, /3. 

Relations (23)-(29) form the full nonlinear MHD problem in the thin disc approxima- 
tion and are named as defined above, the reduced nonlinear equations. 



The linear problem 

Assuming now that the perturbations are small, the system of equations (23)-(28) 
may be linearized about the steady-state equilibrium solution. The resulting system of 
equations is composed of two uncoupled systems, namely which describe the Alfven- 
Coriolis (AC) and the Magnetosonic (MS) modes, respectively. 



Linear stability analysis for the Alfven-Coriolis modes. 

We start by representing the perturbations up to a radius-dependent amplitude factor 
as follows: 

f(r,ri,t) = f(r,ri)cxp[-a(r)d(r)t], (30) 
where X (r) is the complex eigenvalue 

A=A+/r. (31) 



Substituting (30)-(31) into the linearized problem results in the following system 
of linear ordinary differential equations for the perturbed velocity as well as in-plane 



magnetic field components. That system of equations characterizes the Alfven-Coriolis 
waves and depends parametrically on the radius: 



1 db r 

-iAv r - 2v e - 5 — - -7- = 0, (32) 

-iAv e + -v r - - =0, (33) 

2 j3(r)v(Tj) dr\ 

dv r 

-iXb r - = 0, (34) 

-ahe- d -£ + 3 -b r = 0. (35) 
arj 2 

In addition, the Alfven-Coriolis sub-system (32)-(35) is subject to the vanishing bound- 
ary conditions for the in-plane magnetic-field components. 

The linear set of equations (32)-(35) may be reduced to the following single fourth 
order ordinary differential equations for both b r and be: 



d 

dr\ 



1 d 2 I 1 db r ,e 



v{r\) dr\ 2 \ v(t]) dr\ 



H3 + 2l>Mr)± f^^j +lHl 2 -mr)ke = 0. (36) 

Equation (36) is the same as the one used by [5] who have derived it for a model 
problem under the assumption of zero radial variations of the perturbations. Here, it 
is important to emphasize however that the radial coordinate is a parameter a fact that 
renders the radial dependence of the perturbations arbitrary. [5] have solved (36) with the 
aid of the Wentzel-Kramers-Brillouin (WKB) approximation for v(r]) = exp(— f] 2 /2). 
Remarkably, however, a full analytical solution of Eq. (36) is possible for a slightly 
modified density profile. A detailed description of the solution may be found in [9]. We 
repeat here the main results. The first and main step towards that goal is to replace 
the isothermal density vertical steady-state distribution v(r/) = exp(— r/ 2 /2) by the 
following function: 

v(tj) = sech 2 (M]), (37) 

where the shape parameter b is determined by the requirement that the total mass of the 
disc does not change, namely: 

poo poo 

/ exp(-T] 2 /2)d?7 = / sech 2 {br})dt}. (38) 
Jo Jo 



The result is: 



b=^ljn. 



(39) 



As shown in [9] the results derived by employing that profile are hardly distinguishable 
from the WKB results obtained for the true exponential distribution. In fact, notwith- 
standing the use of the terms true and model profiles, such a change of the equilibrium 
profile of the number density (that is determined by the axial momentum balance equa- 
tion) may actually represent some true equilibrium that is obtained from a slightly dif- 
ferent gravitational potential (see [10], where a similar density profile has been obtained 
as an exact solution for flat disc-galaxies whose disc mass content is larger than the mass 
of the central object). 

Introducing a new independent variable % = tanh(Z?T]), such that — 1 < £ < 1, a simple 
equation emerges, which may be cast into the following form: 

(L + K-)(L + K+)v e =0, (40) 

where L is the Legendre operator of second order: 

L= ^ [(1 "^ ] ' * ± = J#+ 2 ^ 2 ±v / 9Ti6n 

Imposing now the zero boundary conditions of be at r) — > °« leads to the requirement 
that the solution of Eq. (40) for vq diverges polynomially at most when r\ — > °«. It is 
concluded therefore that vq is proportional to the Legendre polynomials Pk{%), and the 
eigenvalues X ± are now determined by the dispersion relation 

i^± = ^l[3+2A 2 ±v / 9 + 16A 2 ] =k(k+l). (41) 
2b 

Setting the arbitrary amplitude of vq to unity, and using the linear equations (32) - (35) 
result in the following expressions for the eigenfunctions that are determined up to an 
arbitrary radius dependent amplitude factor: 

^=^r|[(^ ± -i)(^ ± +i)i8-3][fl k -i(«)-fl k+ i(«)], 
%± ^+1) ^ (1-^)^+1^-1 

6 = ~2k + T4 Ot [Pk-iM-Pk+iM], (42) 

where k = 1,2,... plays the role of axial wave number, b^ = for t, = ±1, since 

i\(l) = 1 and Pk(-1) = (-l) fc for all*. 

Turning back to the dispersion relation (41), it may be written as follows: 

(A ± ) 4 i8 2 -(A ± ) 2 J 8( J 8+6)+9(l- J 8)=0, j3 = JL (43) 
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It is thus obvious that the k-th mode is destabilized when the beta value crosses from 
bellow the threshold that is given by: 



= —Kk+\). 

5K 



(44) 



As a result, a universal (for all values of j8(r) and k) criterion for instability emerges 
which reads: 0(r) > 1. Written in terms of the scaled plasma beta the dispersion 
relation (43) has the following solutions for the eigenvalues of the Alfven-Coriolis 
modes (see Fig. 1): 
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0+6± ^(0+6)2-36(1-0) 



2/3 



(45) 



The two eigenvalues, A + and A - , represent fast and slow Alfven-Coriolis waves. While 
the fast Alfven-Coriolis modes are always stable, the number of unstable slow modes is 
determined by the plasma beta. The eigenvalues of the slow Alfven-Coriolis modes, A~ 
, are given therefore by: 



a: 



± 



0+6-^(0+6)2-36(1-0) 
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,Im{A~)=0 forj8<l, 



\ 



(/3+6) 2 + 36(/3-l)-/3 -6 



2j8 



Re{A~) = for > 1. (46) 



The eigenvalues are imaginary, and the system is spectrally unstable if > 1, and 
real for < 1 in which case the system is stable. In particular, the minimal critical 
unsealed plasma beta that is needed for instability is determined by the first unstable 

slow Alfven-Coriolis mode, k = 1, and is given by C ^ = 0.42. The unstable modes 
are the well known MRIs and from Eq. (44) it is easy to calculate how many of them 
are excited for a given value of the plasma beta. Thus, there are k unstable modes for 

$cP < P(r) < 0cr +1 ' > - In particular, Eq. (44) yields approximately the square-root law 
for the number of the unstable modes as a function of the plasma beta: 



k<J3KP(r)/2. 



(47) 



Relation (44) or its simplified version (47) is significant for the consequent modeling 
of non-linear development of the instability. It is finally emphasized that the stability 
criterion as well as the number of unstable modes depend on the radius. Thus, different 
areas within the disc may be characterized by different stability properties as well as 
different number of unstable modes. 

The family of the fast Alfven-Coriolis modes, A + (/3), is characterized by frequencies 
that are much larger than the Keplerian frequency for small values of the scaled plasma 



FIGURE 1. Growth rates ±r„ (dashed lines) for the unstable Alfven-Coriolis (MRI) modes and 
frequencies ±A* (solid lines) for the Alfven-Coriolis oscillations vs universal scaled plasma beta, j3 = 
j3(r)/0, = 2k{k+l)/{3n) for the model number density V = sech^^/^Tf); k = 1,2,... is 
the axial wave number. Meshed straight-line asymptotes at j3 ^> 1 are the scaled Keplerian frequencies, 
±A± = ±1. 



beta (large values of axial wave number or small plasma beta) and tends to the Keplerian 
value at large plasma beta values: 



A+=A+ = ± 



J3+6+ v /( J 8+6) 2 -36(l -jS) 
\ V I — , Im{A+) = 0. (48) 



Expressed in terms of the scaled plasma beta, a single figure depicts all possible 
stable as well as unstable modes. This is shown in Fig. 1. The maximal growth rate 
for the unstable modes is achieved around /3 3, which for a given plasma beta value 
determines the axial wave number of the fastest growing modes. 

An illustration of the perturbed toroidal magnetic field is presented in Fig. 2. The 
perturbations are indeed localized within the effective height of the disc which weakly 
depends on the axial wave number. This corresponds to a finite distance between the 
turning points - the natural characteristics of the problem solution in the WKB approx- 
imation (see [5]). 

In order to compare the current results to well known results for infinite homogeneous 
cylinders [3] the growth rate of the unstable modes is depicted in Fig. 3 as a function of 



(a) 



(b) 



FIGURE 2. The toroidal component of the perturbed magnetic field bg vs self-similar axial variable 
77 = z/H(r) for the Legendre polynomials with the axial wave numbers (a) k = 1,2,3,4 and (b) k = 6, 15. 

All curves are calculated for the fixed value of J3 Z = 1.5 (ft = fi z /0 , = b 2 k(k+l)/3, dashed and 
solid curves correspond to the odd and even k, respectively). 

the axial wave number K, where 



L a = V a / Q. is the Alfven length scale, z = Hy/2 is the effective height of the diffused 
disc at which the equilibrium number density, n ~ exp[— z 2 /(2H 2 )], falls by factor 
e . The effective wave number K is the discrete thin-disc analog of the continuous 
wave number for infinite cylindrical discs. For fixed value of the local plasma beta, 
/3 = 0.41, 0.5, 1.5, 2.5,500, the discrete set of the points in the plane {^,r a } is 
presented by one of the interpolating curves 1,2,3,4,5, respectively. The number of 
the discrete points on each interpolating curves corresponds to the admissible values of 

the Alfven-Coriolis mode number k = 1,2,3,... for which j5cP < < j8 (r). Also, the 
range of unstable ^-values is widening as the value of j8 is increased. In particular, at 
large plasma beta the corresponding set of the points due to their large number (curve 
5) should tend to the continuous curve for infinite cylinder geometry in [3]. For finite 
plasma beta there is a discrete number of points on each curves in Fig. 3, where the left 
bound corresponds to the first Alfven-Coriolis mode, k = 1 . 

Thus, for beta values close to ficr\ the number of unstable modes is small and the 
disk stability properties significantly deviate from those predicted by the infinite cylinder 
model. 



'Cl 



k 



K = k 




(49) 



Linear stability problem for the magnetosonic modes 



The magnetosonic sub-system of equations is analyzed now under the same model 
density profile that was employed in the previous section, namely v(t]) = sech 2 (&7]). 



FIGURE 3. Growth rates T a for the slow Alfven-Coriolis modes vs effective wave number K = -jj*j= = 
—4= , k — 1,2,3,... is the number of the Alfven-Coriolis modes, calculated for the model number density 

V = sech(bri). Interpolating curves 1,2,3,4,5 correspond to j3 = 0.41, 0.71, 1.5, 2.5, 500, respectively. 



This unifies the treatment of all linear modes in the system and sets the base for further 
nonlinear analysis. Thus, assuming the form (31) for the perturbed variable, expressed in 
terms of the new independent variable, ^ = tanh(bri), the magnetosonic system [i.e., the 
lineraized version of eqs. (25) and (26)] may be reduced to the following single ordinary 
differential equation for the perturbed number density: 

(l-^ 2 )0 + [T^2+ 2 ]v = O, (50) 

subject to the boundary conditions v(±l) =0. Transforming to a new dependent vari- 
able, 

m = ^^f^) 7 (5D 

the solution to eq. (50) is given by: 

m = V / l^[C 1 / 1 (^) +C 2 / 2 (£)], (52) 

where, 

h = (^fr /2 (M-^) (53) 



FIGURE 4. Two independent solutions of eq. (50) for A = 2. The full (broken) line corresponds to the 
solution with f\ {fi). 

h = (^|r /2 (M + ^), (54) 

and /i = Vl — A 2 . Solutions that satisfy the boundary conditions exist only for A 2 > 0. 
This means that the magnetosonic modes are stable. Furthermore, the magnetosonic 
spectrum that corresponds to the model density profile is continuous. The two indepen- 
dent solutions for A = 2 are depicted in Fig. 4 where it is easy to see that the arbitrary 
(positive) value of A 2 actually determines the axial wave number. 

WEAKLY NONLINEAR ANALYSIS OF THE MRI 

The viability of the MRI as a generator of turbulence in thin discs, and consequently as 
an important cause of angular momentum transport, depends of course on its nonlinear 
dynamical development and saturation mechanisms. In order to gain insight and to 
provide guidelines to full nonlinear numerical simulations we carry out in this section 
a weakly nonlinear analysis of the interaction of the unstable MRI with the stable 
magnetosonic modes. 

A piece-wise constant model 

The eigenvalues as well as the eigenfunctions of the linearized MHD thin-disk system 
of equations, serve as building blocks of the weakly nonlinear analysis of the unstable 



modes, namely the MRI. As was shown in the previous section the eigenfunctions may 
be expressed in terms of the Legendre as well as some well defined hypergeometric 
functions. However, in order to simplify the vast amount of algebra that usually accom- 
panies any attempt to deviate from linearity, a simpler piece-wise constant model of the 
steady state is introduced. As will be shown, the resulting spectrum is very close to that 
obtained from the diffused profiles (see eq. (37)). 

We start by introducing the following piece- wise constant axial profile for the number 
density: 

*(">={o Vo 1:1;!, <*> 

where Vq = yfii/2. 



The spectrum of the AC modes 



Inserting expression (55) into eqs. (32)-(35), results after some algebra in the follow- 
ing dispersion relation for the AC modes: 

[K 2 (k+ - 2A 2 j8] [7t 2 (k+ l -f - (3 + A 2 )j8] -4A 2 j8 2 = 0, k = 0, 1,2.., (56) 

where /3 = Vq/3, and the eigenfunctions are given by: 



b r {r}),b e {T}) 

Mi) Mi) 



(b r0 ,beo)cos[n(k+-)ri] 
(vro,veo)sin[^(fc+^)r]] 



(57) 
(58) 



As in the diffused case discussed above, also for the piece-wise constant profiles the 
number K of unstable modes (i.e. with T\ = —A 2 > for k < K) depends on the value 
of /3. For K — 1 the threshold for the instability is /3 C = K 2 / 12. The latter corresponds to 
0.65 for the diffused profiles (i.e. after eq. (46)) which is close to the exact value. 



The spectrum of the magnetosonic modes 

Expression (55) is inserted now into the linearized version of eqs. (25) and (26). 
Using representation (30) for the perturbations, and utilizing the fact that the perturbed 
axial velocity is an antisymmetric function about the midplane, results in the following 
dispersion relation for the magnetosonic modes: 

X 2 = K 2 (k+ 1 -) 2 , * = 0,1,2..., (59) 
and the following eigenfunctions: 



v z = a k sm[n(k+-)r[} 



(60) 



v = v a k cos[n{k+lj)ri} . (61) 



Weakly nonlinear analysis 

We consider now a beta value which is slightly above j8 c . In that case there is only 
one MRI mode (k = 0), which is characterized by a small growth rate T. The relation 
between T and /3 may be inferred from the dispersion equation. (56). The result may be 
written in the following way: 

j8 = ^(l + « 2 r 2 ), (62) 

where a 2 = 7/9. Thus, up to a radius-dependent multiplicative factor, the eigenfunctions 
of the first unstable mode (k = 0) are given by: 

b r (ri,t) = a(t)(- + — r 2 ) cos -77 (63) 
Ik 7t 

b e (ri,t) = -a(;)_rcos-T7, (64) 

where a(t) = eW r )' is the initially small amplitude of the MRI. Similar expressions may 
be written for v r {r\,t) and VQ{r\,t). Returning now to eq. (25) and (26), we see that the 
MRI mode nonlinearly excites a magnetosonic wave. To lowest order in T and a, the 
forced part of the latter is described, therefore, by: 

pf = -ja 2 (cos7n7 + 1) (65) 

2 o^/Sin^Ti x 
v zJ = -a 2 T(^^ + r 1 ). (66) 

Obviously the excited magnetosonic waves feed back on the MRI. 

The aim of the weakly nonlinear analysis is to find an ordinary differential equation 
that describes the time evolution of the amplitude a(t) of the single unstable MRI mode. 
In order to do that it is realized first that the transition to instability (when /3 reaches the 
value j5 c from below) occurs when the linearized system has a double zero eigenvalue. 
This situation is known as the Takens-Bogdanov bifurcation [1 1]. As a result the sought 
after equation is expected to be of second order as opposed to first order equations (like 
the Landau-Ginzburg one) that characterize systems that bifurcate through a simple zero 
eigenvalue. Thus following Arter [12] it is conjectured that the typical Takens-Bogdanov 
amplitude equation for ideal MHD is of the form: 

— T = T 2 a-aa i . (67) 
at 1 



In order to calculate a it is noticed that the equilibrium solution of eq. (67) provides an 
asymptotic approximation of the steady state solutions of eqs. (23-28). The latter may be 



obtained by employing the Poincare-Lindstedt method. The latter yields the following 
expression for the steady state amplitude of the radial component of the perturbed 
magnetic field: 



Equating the right hand side of eq. (67) to zero and inserting expression (68) results in 
the following expression for a: 

All is ready now to examine the solutions of eq. (67). As an example consider the case 
r = 0.01, a(0) = 0.001, and d(0) = Ta(0). The results of the numerical solution of eq. 
(67) for that particular case are shown in Fig. 5. For small values of a(t) it behaves 
according to the linearized equation, however as the nonlinear term kicks in it varies 
on a much faster time scale, thus giving it a bursty appearance. The saturation level is 
about one order of magnitude higher than the initial condition. In order to examine the 
viability of the weakly nonlinear analysis the results presented in Fig. 5 are compared 
to the results obtained form the numerical solution of the fully nonlinear set of reduced 
equations (23-28). The initial conditions are given by eqs. (64) and (66) while /3 is 
determined from eq. (62) with Y = 0.01. That solution for the time evolution of the 
amplitude of the radial component of the perturbed magnetic field is displayed in Fig. 
6. It is clear that the weakly nonlinear analysis provides an accurate description of the 
evolution of the MRI near the threshold. Numerical solutions for initial conditions that 
are far from threshold indicate that the nature of oscillatory saturation of the instability 
is preserved. In that case however, the saturation occurs on realistic time scale of the 
order of ten rotation times. 



CONCLUSIONS 

The full spectrum of the MHD modes in thin rotating axially-isothermal discs under 
the influence of axial magnetic field has been obtained analytically. The number of 
unstable MRI modes as a function of the plasma beta has been derived, as well as 
the dispersion relation for the stable magnetosonic waves. A weakly nonlinear analysis 
points out one possible mechanism of saturation of the MRI, namely energy transfer 
from the latter to the stable magnetosonic modes. For the current model this occurs 
as oscillatory saturation due to the conservative nature of the Ideal MHD equations. 
Adding enough dissipation, the degeneracy of the bifurcation point may be removed, 
which ultimately leads to a first order Landau-Ginzburg like equation for the amplitude, 
instead of the second order equation given in eq. (67). Saturation of the latter type has 
been investigated by Umurhan et al. [13]. 
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